Principal Component Analisis

Análisis de Componentes Principales con DESeq2

Instalar y cargar la librería DESeq2. Esta librería contiene una serie de funcionalidades para realizar análisis estadísticos, “aguas abajo” (down-stream) de datos ómicos analizados con SALMON e importados con tximport.

Importar las librerias necesarias

Son necesarias las librerías de “DESeq2”, “readxl”, “ggplot2”, “ggrepel” y “plotly”.

# Instalar, si no está instalado, BiocManager: cargador de dependencias de bioconductor
if (!require("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

## Warning: package 'BiocManager' was built under R version 4.2.1

# Instalar, si no está instalado, DESeq2: herramienta para el análisis transcriptómico aguas abajo
if (!require("DESeq2", quietly = TRUE)) {
  BiocManager::install("DESeq2")
}

## Warning: package 'DESeq2' was built under R version 4.2.2

## Warning: package 'S4Vectors' was built under R version 4.2.1

## Warning: package 'BiocGenerics' was built under R version 4.2.1

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Warning: package 'IRanges' was built under R version 4.2.1

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:grDevices':
## 
##     windows

## Warning: package 'GenomeInfoDb' was built under R version 4.2.2

## Warning: package 'SummarizedExperiment' was built under R version 4.2.1

## Warning: package 'MatrixGenerics' was built under R version 4.2.1

## Warning: package 'matrixStats' was built under R version 4.2.2

## 
## Attaching package: 'MatrixGenerics'

## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars

## Warning: package 'Biobase' was built under R version 4.2.1

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians

## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

# Instalar, si no está instalado, weitexl: paquete para guardar como excel dataframes
if (!require("writexl", quietly = TRUE)) {
  install.packages("writexl")
}

## Warning: package 'writexl' was built under R version 4.2.2

# NOTA: si dice de actualizar paquetes, contestar que NO (n).

# Instalar, si no está instalado
# - readxl: sirve para leer excele
# - ggplot2: sirve para generar gráficos estáticos, muy usado en publicaciones
# - ggrepel: sirve para añadir funciones a ggplot2 que evitan que los textos dentro de un gráfico se superpongan
# - plotly: sirve para generar gráficos dinámicos, permite interactuar al usuario con la gráfica
# Normalmente, vienen instalado en muchas distribuciones de R, solo hace falta cargarlos con "library()"

Leer los datos

Cargamos los datos clínicos de los pacientes y los datos del transcriptoma.

# --- Datos de expresion genética ---

# El objeto que contiene toda la información esta guardado como .RDS
txi <- readRDS("../import quant/txi.RDS")
# Trabajaremos con los counts
countData <- txi$counts

# No vamos a trabajar más con txi, así que como ocupa 35MB, lo elimino de la memoria
rm(txi)

# La matriz tiene, en los nombres de las muetras, un sufijo después de una barra baja, que es el "sample ID" que se utilizó para secuenciarlo en illumina. Para que no de problemas, se van a renombrar

# Elimino el sufijo "_Sxx" de los nombres de las columnas
lNewColnames <- unlist(strsplit(colnames(countData), "_"))[seq(1, 2 * ncol(countData), 2)]

# Asigno los nuevos nombres de las columnas a la matriz
colnames(countData) <- lNewColnames

# Muestro los primeros datos
head(countData)

##                        103     105     112     113     131     151     165
## ENSG00000000003.15 280.690 269.872 437.320 354.990 565.889 474.571 393.377
## ENSG00000000005.6    0.000   1.000   2.000   5.000  14.000   0.000   1.000
## ENSG00000000419.14 468.499 657.489 417.199 670.997 476.536 303.273 464.651
## ENSG00000000457.14 292.789 289.795 250.316 338.842 339.367 258.275 297.864
## ENSG00000000460.17 148.557 101.415 201.240 132.638 132.400 169.704 270.774
## ENSG00000000938.13  37.000  84.008  82.001  36.000  87.000 108.000  39.180
##                         17     171     177     188     201     203     210
## ENSG00000000003.15 537.758 248.298 295.417 678.890 518.671 377.485 602.191
## ENSG00000000005.6    4.000   1.000   0.000   0.000   0.000   0.000   1.000
## ENSG00000000419.14 353.991 540.465 459.604 386.604 494.334 247.871 623.557
## ENSG00000000457.14 259.594 223.485 276.667 319.040 335.244 197.801 352.266
## ENSG00000000460.17 159.832 141.252 192.590 143.592 247.313 121.909 272.588
## ENSG00000000938.13  82.001  37.000  49.000  38.000 103.999 114.000  57.000
##                        218      30      37       5      63      64      70
## ENSG00000000003.15 636.489 418.424 410.440 597.728 460.811 463.123 625.235
## ENSG00000000005.6    2.000   3.000   0.000   0.000   1.000   1.000   1.000
## ENSG00000000419.14 616.128 881.789 460.925 457.024 566.649 324.236 852.521
## ENSG00000000457.14 370.288 464.448 238.741 257.118 287.889 256.463 359.007
## ENSG00000000460.17 302.957 212.938 222.348 178.380 145.989 125.019 226.644
## ENSG00000000938.13  18.000  88.000  70.999 144.001  44.000  59.000 144.590
##                         77      78      80      83       88      95      98
## ENSG00000000003.15 438.495 425.433 250.667 416.714 3950.941 661.838 572.685
## ENSG00000000005.6    0.000   2.000   4.000   0.000   15.000   0.000   0.000
## ENSG00000000419.14 605.586 323.024 824.963 543.388 4038.405 825.350 526.198
## ENSG00000000457.14 313.623 232.061 315.194 242.842 2167.094 286.643 397.920
## ENSG00000000460.17 207.402 214.585 109.392 204.298 1210.723 161.989 267.191
## ENSG00000000938.13  69.000  91.030 244.000  47.001  277.017  92.343  77.000

rm(lNewColnames)

# --- Datos clínicos ---

# Cargar la librería necesaria para leer archivos excel
library("readxl")

## Warning: package 'readxl' was built under R version 4.2.1

# Leer los datos del excel
clinicalData <- read_excel("../datos clinicos/dfClinical.xlsx")

# Muestro los primeros datos
head(clinicalData)

## # A tibble: 6 × 21
##   Samples Age.at.diagnosis Gender Tumor.localization  Survival Grade Ki67  LOH  
##     <dbl>            <dbl> <chr>  <chr>               <chr>    <chr> <chr> <chr>
## 1     177               63 F      Left frontal lobe   118      IV    N.D.  0    
## 2      80               71 F      Left frontal lobe   244      IV    25    0    
## 3      88               45 M      Left frontal lobe   3384     II    5     N.D. 
## 4       5               53 M      Left frontal lobe   1239     IV    N.D.  0    
## 5      70               67 M      Left fronto-pariet… 406      IV    N.D.  0    
## 6      83               52 F      Left parietal lobe  22       IV    80    N.D. 
## # … with 13 more variables: `IDH1.(DM)` <chr>, `IDH2.(DM)` <chr>, MGMT <chr>,
## #   pTERT.C228T <chr>, pTERT.C250T <chr>, NLR <chr>, PLR <chr>, EGFRvIII <chr>,
## #   fAge <chr>, fGender <chr>, fTumorHemisphere <chr>, fTumorLobule <chr>,
## #   fSurvival <chr>

Para evitar errores a la hora de fusionar los datos de la cuantificación y los datos de los pacientes, se va a reordenar el dataframe de los datos clínicos de los pacientes de forma que tengan el mismo orden que las columnas de la matriz con los counts para cada gen.

# Los datos de la matriz estan ordenados como si los números fuesen carácteres, por lo que voy a transformar la variable "Sample" del dataframe de los datos clínicos a carácter, y después lo ordenaré de menor a mayor
clinicalData$Samples <- as.character(clinicalData$Samples)
clinicalData <- clinicalData[order(clinicalData$Samples),]

Elimino los datos del paciente 88, ya que por un error en el laboratorio, se secuenció con más profundidad de lectura y sus datos de supervivencia son outliers.

# Elimino la columna del paciente 88 de la matriz "countData" 
countData <- countData[,colnames(countData) != "88"]

# Elimino la fila del paciente 88 en los datos clinicos
clinicalData <- clinicalData[clinicalData$Samples != "88",]

# Nos aseguramos de que hay el número correcto de muestras
ncol(countData)

## [1] 27

nrow(clinicalData)

## [1] 27

# Mostramos
head(countData)

##                        103     105     112     113     131     151     165
## ENSG00000000003.15 280.690 269.872 437.320 354.990 565.889 474.571 393.377
## ENSG00000000005.6    0.000   1.000   2.000   5.000  14.000   0.000   1.000
## ENSG00000000419.14 468.499 657.489 417.199 670.997 476.536 303.273 464.651
## ENSG00000000457.14 292.789 289.795 250.316 338.842 339.367 258.275 297.864
## ENSG00000000460.17 148.557 101.415 201.240 132.638 132.400 169.704 270.774
## ENSG00000000938.13  37.000  84.008  82.001  36.000  87.000 108.000  39.180
##                         17     171     177     188     201     203     210
## ENSG00000000003.15 537.758 248.298 295.417 678.890 518.671 377.485 602.191
## ENSG00000000005.6    4.000   1.000   0.000   0.000   0.000   0.000   1.000
## ENSG00000000419.14 353.991 540.465 459.604 386.604 494.334 247.871 623.557
## ENSG00000000457.14 259.594 223.485 276.667 319.040 335.244 197.801 352.266
## ENSG00000000460.17 159.832 141.252 192.590 143.592 247.313 121.909 272.588
## ENSG00000000938.13  82.001  37.000  49.000  38.000 103.999 114.000  57.000
##                        218      30      37       5      63      64      70
## ENSG00000000003.15 636.489 418.424 410.440 597.728 460.811 463.123 625.235
## ENSG00000000005.6    2.000   3.000   0.000   0.000   1.000   1.000   1.000
## ENSG00000000419.14 616.128 881.789 460.925 457.024 566.649 324.236 852.521
## ENSG00000000457.14 370.288 464.448 238.741 257.118 287.889 256.463 359.007
## ENSG00000000460.17 302.957 212.938 222.348 178.380 145.989 125.019 226.644
## ENSG00000000938.13  18.000  88.000  70.999 144.001  44.000  59.000 144.590
##                         77      78      80      83      95      98
## ENSG00000000003.15 438.495 425.433 250.667 416.714 661.838 572.685
## ENSG00000000005.6    0.000   2.000   4.000   0.000   0.000   0.000
## ENSG00000000419.14 605.586 323.024 824.963 543.388 825.350 526.198
## ENSG00000000457.14 313.623 232.061 315.194 242.842 286.643 397.920
## ENSG00000000460.17 207.402 214.585 109.392 204.298 161.989 267.191
## ENSG00000000938.13  69.000  91.030 244.000  47.001  92.343  77.000

head(clinicalData)

## # A tibble: 6 × 21
##   Samples Age.at.diagnosis Gender Tumor.localization  Survival Grade Ki67  LOH  
##   <chr>              <dbl> <chr>  <chr>               <chr>    <chr> <chr> <chr>
## 1 103                   63 F      Right fronto-parie… 131      IV    N.D.  0    
## 2 105                   62 M      Left temporal lobe  648      IV    N.D.  0    
## 3 112                   53 M      Right occipital lo… 679      IV    N.D.  N.D. 
## 4 113                   45 M      Right temporal lobe N.D.     II    N.D.  0    
## 5 131                   68 M      Right temporal lobe 23       IV    60    N.D. 
## 6 151                   46 F      Right parieto-temp… 73       IV    N.D.  N.D. 
## # … with 13 more variables: `IDH1.(DM)` <chr>, `IDH2.(DM)` <chr>, MGMT <chr>,
## #   pTERT.C228T <chr>, pTERT.C250T <chr>, NLR <chr>, PLR <chr>, EGFRvIII <chr>,
## #   fAge <chr>, fGender <chr>, fTumorHemisphere <chr>, fTumorLobule <chr>,
## #   fSurvival <chr>

Aplicar la transformación de DESeq2

La librería DESeq2, permite fusionar los datos de los transcritos y los datos clínicos de los pacientes. De forma que por un lado tenemos la matriz con los datos de expresion génica (counts de los transcritos) y por otro lado tenemos los datos clínicos de los pacientes que son tratados como fáctores.

@see: “https://compbiocore.github.io/deseq-workshop-1/assets/deseq_workshop_1.html”

# Cargar la librería
library("DESeq2")

# Aplicar la transformación que une counts + metadados
dds <- DESeqDataSetFromMatrix(countData = round(countData), colData = clinicalData, design = ~Grade)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

# Mostramos la información del objeto (es un objeto especial de DESeq2)
dds

## class: DESeqDataSet 
## dim: 38244 27 
## metadata(1): version
## assays(1): counts
## rownames(38244): ENSG00000000003.15 ENSG00000000005.6 ...
##   ENSG00000291240.1 ENSG00000291292.1
## rowData names(0):
## colnames(27): 103 105 ... 95 98
## colData names(21): Samples Age.at.diagnosis ... fTumorLobule fSurvival

PCA con DESeq2

Ahora vamos a aplicar un Analisis de Componentes Principales (PCA) con las funciones del paquete DESeq2.

# Normalizamos los datos
# vst = Variance Stabilizing Transformation
vsdata <- vst(dds, blind = FALSE)

# Realizamos la extracción de componentes principales
pcaData <- plotPCA(vsdata, intgroup="Grade", returnData = TRUE)

# plotPCA, por defecto, genera un dataframe con una variable llamada "group" que contiene la misma información que la variable que se ha indicado en la función. Se pude eliminar.
pcaData$group <- NULL

El objeto “pcaData” contiene los valores de las 2 primeras Componentes Principales (PCs) para los 28 pacientes.

Se ha creado una función, partiendo de la original, que calcula las 4 primeras componentes principales.

# Mostramos los datos
head(pcaData)

##           PC1         PC2 Grade name
## 103 -46.78127   5.9369410    IV  103
## 105 -50.69042   6.1853577    IV  105
## 112  12.29899 -12.6832694    IV  112
## 113 -55.66857  -0.6387503    II  113
## 131  10.82648  14.7033023    IV  131
## 151  16.34288  20.5112624    IV  151

Grafico de dispersión

# La ubicación donde guardar el gráfico
dir.create(file.path(getwd(), "img"), showWarnings = FALSE)

# Cargamos la librería ggplot2
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.2

library(ggrepel)

## Warning: package 'ggrepel' was built under R version 4.2.1

# Extraemos el atributo del porcentaje de la varianza
# (sirve para añadir el % de la varianza en los ejes del gráfico)
percentVar <- round(100 * attr(pcaData, "percentVar"))

# Hacemos el gráfico
g <- ggplot(pcaData, aes(PC1, PC2, color=Grade, label = row.names(pcaData))) +
  geom_point(size=3) +
  geom_label_repel(size = 3) + 
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  ggtitle("Patients by glioblastoma grade") + 
  coord_fixed();g

# Guardamos el gráfico
ggsave("./img/PCA_01_fGrade.png", plot = g)

## Saving 7 x 5 in image

rm(percentVar, g)

PCA CP3 y 4

La función “plotPCA” del paquete DESeq2 sola genera las Componentes Principales (CP) 1 y 2. Como Luis Miguel represento los datos de los pacientes con 3 CP, es necesario calcular las otras componentes principales.

# ----------- #
# --- PCA --- #
# ----------- #

# --- Variables --- #
# Número de genes
ntop = 500 
# El fáctor
intgroup = "Grade"

# --- Extracción de PCs --- #
rv <- rowVars(assay(vsdata))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca <- prcomp(t(assay(vsdata)[select, ]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)

intgroup.df <- as.data.frame(colData(vsdata)[, intgroup, drop = FALSE])

# --- Dataframe con las PCs --- #
pcaData2 <- data.frame(
  PC1 = pca$x[, 1], 
  PC2 = pca$x[, 2], 
  PC3 = pca$x[, 3],
  PC4 = pca$x[, 4], 
  intgroup.df, 
  name = colData(dds)[,1])

# Mostramos el % de la varianza original retenida (con la varianza se estima la información)
round(sum(percentVar[1:4]) * 100, 2)

## [1] 57.31

# Eliminar variables
rm(ntop, intgroup, rv, select, pca,  intgroup.df)

# Mostramos el dataframe generando con las componentes principales 1 - 4
head(pcaData2)

##           PC1         PC2         PC3        PC4 Grade name
## 103 -46.78127   5.9369410 -14.2498829   8.051700    IV  103
## 105 -50.69042   6.1853577  -6.6862629  -1.714114    IV  105
## 112  12.29899 -12.6832694  -0.1276191  -6.318918    IV  112
## 113 -55.66857  -0.6387503  -6.7935896  -2.988240    II  113
## 131  10.82648  14.7033023  13.7745312  16.942380    IV  131
## 151  16.34288  20.5112624   6.2481519 -15.658440    IV  151

PCA con más variables

Por defecto, la función que realiza el PCA del paquete DESeq2 selecciona los primeros 500 genes en función de su varianza. Es decir, ordena de mayor a menor varianza las variables y selecciona las primeras 500. La poca retención de información de las componentes principales (CP) se debe en parte a que se están tomando una muestra pequeña de las variables originales (de alrededor de 37.000, solo se trabaja con 500) y porque la relación entre las variables posiblemente no sea líneal.

# ----------- #
# --- PCA --- #
# ----------- #

# --- Variables --- #
# Número de genes
ntop = 2000 
# El fáctor
intgroup = "Grade"

# --- Extracción de PCs --- #
rv <- rowVars(assay(vsdata))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca <- prcomp(t(assay(vsdata)[select, ]))
percentVar3 <- pca$sdev^2/sum(pca$sdev^2)

intgroup.df <- as.data.frame(colData(vsdata)[, intgroup, drop = FALSE])

# --- Dataframe con las PCs --- #
pcaData3 <- data.frame(
  PC1 = pca$x[, 1], 
  PC2 = pca$x[, 2], 
  PC3 = pca$x[, 3],
  PC4 = pca$x[, 4], 
  intgroup.df, 
  name = colData(dds)[,1])

# Mostramos el % de la varianza original retenida (con la varianza se estima la información)
round(sum(percentVar3[1:4]) * 100, 2)

## [1] 57

# Eliminar variables
rm(ntop, intgroup, rv, select, pca,  intgroup.df)

# Mostramos el dataframe generando con las componentes principales 1 - 4
head(pcaData3)

##           PC1       PC2          PC3        PC4 Grade name
## 103 -69.46740 -10.39328   2.69085208 -21.063693    IV  103
## 105 -75.02701 -16.21282  -0.07726929 -17.733730    IV  105
## 112  17.23165  11.66434  -1.66402445   6.250943    IV  112
## 113 -84.04938  -3.62071  -4.26842827 -15.594621    II  113
## 131  14.63952 -14.14128  37.05831890  12.940433    IV  131
## 151  26.65543 -37.64575 -20.23258993  20.838870    IV  151

# Guardamos los datos
# write_csv(pcaData3, "./PCA_03.csv")

# Cargamos la librería ggplot2
library(ggplot2)
library(ggrepel)

# Hacemos el gráfico
g <- ggplot(pcaData3, aes(PC1, PC2, color=Grade, label = row.names(pcaData3))) +
  geom_point(size=3) +
  geom_label_repel(size = 3) + 
  xlab(paste0("PC1: ",round(percentVar3[1] * 100, 2),"% variance")) +
  ylab(paste0("PC2: ",round(percentVar3[2] * 100, 2),"% variance")) + 
  ggtitle("Patients by glioblastoma grade") + 
  coord_fixed();g

# Guardamos el gráfico
ggsave("./img/PCA_03_fGrade.png", plot = g)

## Saving 7 x 5 in image

rm(percentVar, g)

Gráfico de dipersión en 3D

Se utilizará un paquete distinto a ggplot2 para realizar gráficas: plotly. Plotyl es un paquete utilizado para generar gráficos dinámicos, los cuales permiten al usuario interactuar con ellos. Para la visualización de estos datos, es una mejor opción que ggplot2.

# Representación de los individuos con las 3 primeras componentes principales.
# Cargamos la librería plotly
library(plotly)

## Warning: package 'plotly' was built under R version 4.2.1

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:IRanges':
## 
##     slice

## The following object is masked from 'package:S4Vectors':
## 
##     rename

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# Generamos la figura
fig <- plot_ly(
  pcaData2,
  x = ~PC1,
  y = ~PC2,
  z = ~PC3
) %>% add_markers(color = ~Grade); fig

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

# Generamos la figura, con la opción de que cada punto tenga indicado el nombre de la muestra.
fig_named <- fig %>% add_trace(text = row.names(pcaData2), hoverinfo = 'text', marker = list(color='green')); fig_named

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

rm(fig, fig_named)

Añadir los factores de estudio al dataframe con las Componentes Principales

Al dataframe que se ha generado con la función del paquete DESeq2, se le van a añadir los factores que se encuentran en el dataframe de los datos clínicos de los pacientes.

# --- PCA Data 1 --- #
# Factor edad
pcaData$fAge <- clinicalData$fAge
# Factor género
pcaData$fGender <- clinicalData$fGender
# Factor hemisferio donde aparecio el tumor
pcaData$fTumorHemisphere <- clinicalData$fTumorHemisphere
# Factor lóbulo donde aparecio el tumor
pcaData$fTumorLobule <- clinicalData$fTumorLobule
# Factor supervivencia
pcaData$fSurvival <- clinicalData$fSurvival

# --- PCA Data 2 --- #
# Factor edad
pcaData2$fAge <- clinicalData$fAge
# Factor género
pcaData2$fGender <- clinicalData$fGender
# Factor hemisferio donde aparecio el tumor
pcaData2$fTumorHemisphere <- clinicalData$fTumorHemisphere
# Factor lóbulo donde aparecio el tumor
pcaData2$fTumorLobule <- clinicalData$fTumorLobule
# Factor supervivencia
pcaData2$fSurvival <- clinicalData$fSurvival

Gráficos de dispersion según factores

Utilizando los factores de los datos clínicos de los pacientes, podemos comprobar si alguno de ellos muestra un patrón.

Gráficos de tres componentes principales.

Fáctor hemispherio.

# Representación de los individuos con las 3 primeras componentes principales.
# Cargamos la librería plotly
library(plotly)

# Generamos la figura
fig <- plot_ly(
  pcaData2,
  x = ~PC1,
  y = ~PC2,
  z = ~PC3
) %>% add_markers(color = ~fTumorHemisphere); fig

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

rm(fig)

Fáctor lóbulo cerebral.

# Representación de los individuos con las 3 primeras componentes principales.
# Cargamos la librería plotly
library(plotly)

# Generamos la figura
fig <- plot_ly(
  pcaData2,
  x = ~PC1,
  y = ~PC2,
  z = ~PC3
) %>% add_markers(color = ~fTumorLobule); fig

rm(fig)

Gráficos de dos componentes principales.

# Extraemos el atributo del porcentaje de la varianza
# (sirve para añadir el % de la varianza en los ejes del gráfico)
percentVar <- round(100 * attr(pcaData, "percentVar"))

# Xlab y Ylab
sXlab <- paste0("PC1: ", percentVar[1], "% variace")
sYlab <- paste0("PC2: ", percentVar[2], "% variace")

# Factor edad
g1 <- ggplot(pcaData, aes(PC1, PC2, color=fAge, label = row.names(pcaData))) +
  geom_point(size=3) +
  geom_label_repel(size = 3) + 
  xlab(sXlab) +
  ylab(sYlab) + 
  ggtitle("Patients by age quartils") +
  coord_fixed()

# Factor hemisferio cerebral
g2 <- ggplot(pcaData, aes(PC1, PC2, color=fTumorHemisphere, label = row.names(pcaData))) +
  geom_point(size=3) +
  geom_label_repel(size = 3) + 
  xlab(sXlab) +
  ylab(sYlab) +
  ggtitle("Patients by tumor hemisphere") +
  coord_fixed()

# Factor lobulo cerebral
g3 <- ggplot(pcaData, aes(PC1, PC2, color=fTumorLobule, label = row.names(pcaData))) +
  geom_point(size=3) +
  geom_label_repel(size = 3) + 
  xlab(sXlab) +
  ylab(sYlab) + 
  ggtitle("Patients by tumor lobule") +
  coord_fixed()

# Factor lobulo + hemisferio cerebral
g4 <- ggplot(pcaData, aes(PC1, PC2, color=fTumorLobule, shape = fTumorHemisphere, label = row.names(pcaData))) +
  geom_point(size=3) +
  geom_label_repel(size = 3) + 
  xlab(sXlab) +
  ylab(sYlab) + 
  ggtitle("Patients by tumor location") +
  coord_fixed()

# Factor superviviencia
g5 <- ggplot(pcaData, aes(PC1, PC2, color=as.factor(fSurvival), label = row.names(pcaData))) +
  geom_point(size=3) +
  geom_label_repel(size = 3) + 
  xlab(sXlab) +
  ylab(sYlab) + 
  ggtitle("Patients by survival") +
  coord_fixed()

# Mostramos todos los gráficos
g1;g2;g3;g4;g5

# Guardamos los gráficos
ggsave("./img/PCA_02_fAge.png", plot = g1)

## Saving 7 x 5 in image

ggsave("./img/PCA_03_fHemisphere.png", plot = g2)

## Saving 7 x 5 in image

ggsave("./img/PCA_04_fLobule.png", plot = g3)

## Saving 7 x 5 in image

ggsave("./img/PCA_05_fHemisphere_fLobule.png", plot = g4)

## Saving 7 x 5 in image

ggsave("./img/PCA_06_fSurvival.png", plot = g5)

## Saving 7 x 5 in image

rm(percentVar, sXlab, sYlab, g1, g2, g3, g4, g5)

Guardar las Componentes Principales

Guardar como los dataframe con las componentes principales como exceles. También, se guarda como un archivo “.RDS”.

# --- Guardar como Excel --- #
library("writexl")

# Guardamos las componentes principales 1
write_xlsx(pcaData, "./pcaData.xlsx")

# Guardamos las componentes principales 2
write_xlsx(pcaData2, "./pcaData2.xlsx")

# --- Guardar como RDS --- #
saveRDS(pcaData, "./pcaData.rds")
saveRDS(pcaData2, "./pcaData2.rds")

Borramos todas las variables

rm(pcaData, pcaData2, dds, countData, vsdata, clinicalData)